Draft version July 24, 2012 

Preprint typeset using I^'T^]X style emulatcapj v. 8/13/10 



TURBULENCE INDUCED COLLISIONAL VELOCITIES AND DENSITY ENHANCEMENTS: LARGE 
INERTIAL RANGE RESULTS FROM SHELL MODELS 

Alexander Hubbard^'^ 

^ Max Planck Institut fiir Astronomic, Konigstuhl 17, D-69117 Heidelberg, Germany 
2 Department of Astrophysics, American Museum of Natural History, 79th St. at Central Park West, New York, NY 10024-5192, USA 



(N 
O 



m 
(N 

o 

in: 



> 

in 
m 

o 



Draft version July 24, 2012 

ABSTRACT 

To understand the earliest stages of planet formation, it is crucial to be able to predict the rate 
and the outcome of dust grains collisions, be it sticking and growth, bouncing, or fragmentation. The 
outcome of such collisions depends on the collision speed, so we need a solid understanding of the 
rate and velocity distribution of turbulence-induced dust grain collisions. The rate of the collisions 
depends both on the speed of the collisions and the degree of clustering experienced by the dust 
grains, which is a known outcome of turbulence. We evolve the motion of dust grains in simulated 
turbulence, an approach that allows a large turbulent inertial range making it possible to investigate 
the effect of turbulence on meso-scale grains (millimeter and centimeter). We find three populations 
of dust grains: one highly clustered, cold and collisionless; one warm; and the third "hot" . Our results 
can be fit by a simple formula, and predict both significantly slower typical collisional velocities for a 
given turbulent strength than previously considered, and modest effective clustering of the collisional 
populations, easing difficulties associated with bouncing and fragmentation barriers to dust grain 
growth. Nonetheless, the rate of high velocity collisions falls off merely exponentially with relative 
velocity so some mid- or high- velocity collisions will still occur, promising some fragmentation. 

Subject headings: turbulence - planetary systems: protoplanetary disks - planets and satellites: for- 
mation 



1. INTRODUCTION 



Collisions between 
tary disks is a key 
mation, and so has 

of study, bo th theo reti cal (iVolk et al.l 



dust grains in protoplane- 
process in planetesimal for- 
significant amount 



seen a 



Markiewicz et al.l 



119911: 

Youdin fc Goodman 2005; Ormcl fc Cuzzil I2007D and 
nume rical fjPullemond fc Dominiki2005l : IJohansen et al.l 
l2007f l. More, the results of collisions between arti- 
ficial dust grains as a function of relative velocity 
can be directly obse r ved in labora t ory e xperiments 
(|Blum fc Wurmi 12008) : IGuttler et all l20Tol) . These 
collisions can lead to sticking and growth of dust grains, 
or to fragmentation, repopulating the smallest sizes of 
dust grains. On the other hand, laboratory experiments 
suggest that the null-result of a collision, bouncing, 
poses a significant threat to dust growth beyond the cm 
scale (jZsom et al.|[2010, ). 

The gas component of protoplanetary disks is believed 
to be turbulent, and the effect of this turbulence on the 
collisions between dust particles entrained in the flow can 
be reasonably estimated analytica lly, for ex ample as done 
in the line of work starting wit h IVolk et al. ( 19 8Q) and 
more recently elaborated on bv lOrmel fc CuzzI (j2007[ ). 
Such analyses give the unsurprising result that dust col- 
lisional velocities are comparable to the turbulent veloc- 
ities of the gas on scales set by the properties of the 
dust grains. These estimated velocities are, however, ex- 
pected to be large. The turbulence in protoplanetary 
disks is often invoked to explain th e disk's accretion on 
the mega- year time-scales observed (jShakura fc SunvaevI 
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iCuz zi fc Hoga,nl 



1980; 



2003 



Il973t IRussell et al.l 120061 ). and the turbulent motions 
required to achieve such a feat are expected to drive 
dust collisions th at destroy the participants or merely re- 
sult in bou ncing (jWurm et al.ll200lllGuttler et"alll2010l: 
iWettlaufeil |2010[ ). The existence of the bouncing and 
fragmentation velocity cut-offs means that it is important 
to determine not merely the rate of collisions, but also 
the collisional probability distribution as a function of 
the relative velocity, because outlying e vents could cross 
the bouncing barrier (jZsom et al.ll2010[ ). 

Another important behavior of particles in tur- 
bulence 
19871 



is "prefer ent ial con centration" (IMaxeyl 

IFessler et al.l p9l ICuzzi et al.l I200li: 

Toschi fc Bodenschatd '2009), where inertial parti- 



cles are ejected from regions of high vorticity due to 
centrifugal forces, and accumulate in regions of high 
strain. This has been hypothesized to increase dust 
collisi on rates by crea ting local dust density fluctua- 
tions (IPan et aU 120 111 ), as well as possibly leading to 
dust drag on the turbulence. The latter can cause a 
streaming-inst ability that further enhances the local 
dust density (, Youdin fc GoodmanI 120051 : IJohansen et al.l 
120071 ). Beyond the effect on dust grain collisions rates, 
such density enhancements can provide a rapid route to 
gravitational i nstabilities and t he resulting collapse into 
planetesimals (j Johansen et al.l 120061) . This ejection of 
inertial particles from turbulent (intermittent) vorticies, 
which will be discussed in this paper, should be distin- 
guished from the dust-trapping property of anti-cyclonic 
long-lived vor tices, which are large enough to feel 
Corio lis forces (jBarge fc Sommer"ialll995l : IJohansen et al.l 
I2OOI . 

Moreover, the inertial range of the turbulence in these 
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systems is expected to be quite broad, easily five to six 
orders of magnitude in k-space or four orders of magni- 
tude in turbulent turn-over time, as their fluid Reynolds 
numbers are large. This means that there will be dust 
grains well embedded in the turbulent cascade, too small 
to notice that they are trapped within huge eddies, yet 
too large to be meaningfully affected by the smallest 
scale turbulence. These well embedded particles are pre- 
cisely the most important ones, since they include the 
millimeter-centimeter sized dust grains that populate the 
fragmentation and bouncing barriers. A resolution on 
the order of several thousand cubed would be required 
to fully simulate a system evolving the full Navier-Stokes 
equations with an upper double digit inertial range, be- 
yond currently available resources. Further, as we will 
show, turbulent clustering implies that very small par- 
ticle separations must be considered to extract useful 
data. The resolution required for performing a full hy- 
drodynamical simulation is, again, on the order of several 
t housand cubed (S e ction |4T|) . 

iCarballido et al.l ()201Clf ) explicitly run into the above- 
mentioned problems of limited inertial range and low spa- 
tial resolution. We resolve these problems by putting 
particles into a set of turbulent c ascade models known 
as shell models (|Bohr et al.lll998l chapter 3), which al- 
lows us determine the clustering and coUisional velocity 
probability distributions that can be expected for large 
inertial ranges (up to 256) while resolving small particle 
separations. This approach of using synthetic turbulence 
was also used by Bee et al. (2005) , although the details 
of the projection into real space differs. Our work also 
differs by our focus on the collisional velocities as a key 
diagnostic, and we find noticeably lower col hsional ve- 
locities than those estimated analytically bv iVolk et al.l 
(jl980) and following. Further, unlike those works, our 
results are not well approximated by a single effective 
collisional velocity, but require consideration of a veloc- 
ity probability distrib ution. While our clustering results 
are similar to those of iPan et al.l (|2011[ ). simultaneously 
considering the collisional velocities allows us to distin- 
guish between multiple particle collisional populations, 
which changes the conclusions about the effects of the 
strong clustering that is observed. 

In Section [5] we define our numerical methods and in 
Section [3] we use them to investigate a particular base 
case to find the physical behavior. In Section |4] we in- 
vestigate systems with different inertial ranges and syn- 
thetic turbulence implementations. We compare our re- 
sults with previous work in Section [S] and conclude in 
Section [51 We list our main variables, parameters, and 
diagnostics in Table |4] and discuss the definition and use 
of the Stokes number in Appendix 

2. NUMERICAL SETUP 

We use the Pencil CodeQ to track the motion of par- 
ticles in synthetic turbulent cascades. Our synthetic tur- 
bulence derives from shell models: we consider spatial 
Fourier transforms of the turbulent velocity field (from 
X to fe), keeping the modes that fit in a limited range 
of |fe| which defines our considered subsection of the in- 
ertial range. We then coarse-grain the sphere of k vec- 
tors into logarithmically spaced shells in fc-space. As 

^ http:/ /pencil-code. googlecode.com 



such, we consider only the energy spectrum of the tur- 
bulence. We will use four different models . Three fol- 
low a n exact Kolmogorov power spectrum ([Kolmogorovl 
I1941L see Section [2TT|) and differ only in the implementa- 
tion of their time variability, if any. The fourth model, 
discussed in Section [2?TIll combines the spatial approach 
of Section 12.11 with an energy spectrum derived from a 
GOY model, which solves a set of equations for the tur- 
bulent energy cascade. However, it should be noted that 
the magnitude of the velocity field in a GOY model is 
strongly fluctuating in time, making the choice of di- 
mensionless quantities less obvious. A first study with a 
time-independent energy spectrum is needed to establish 
the form of turbulent collisional velocities and is closer 
to existing analytical work. 

Since we do not evolve the Navier-Stokes equations to 
determine the gas motion, we do not need a grid or a sub- 
grid model: this allows for greater resolution. We evolve 
4 X 10^ particles in a periodic box of size (271")^. Our 
setup does not allow for the study of dust backreaction 
on turbulence (preventing, for example, the streaming 
instability), but does allow us to examine dust pairs at 
very small relative separation. 

2.1. Flow field 

Each coarse-grained shell is associated with a fluid ve- 
locity and turnover time and indexed with m which runs 

from TTimin tO TTlmax with 

km = 2™, (1) 
'^m — ^/^mkm- (2) 

In a Kolmogorov spectrum, the velocity obeys the power- 
wave 

Vm = wofc,7i^^^; (3) 

this generates the dashed/red velocity spectrum in Fig- 
ure [1] by construction. Our choice of rrimin and TOmax 
is limited by available numerical resources and by the 
constraints of the periodicity of the numerical box. We 
will vary them for different simulations as the choice de- 
termines the effective inertial range of our turbulence- 
simulate. This allows us to find the effective inertial 
range large enough that the particle response becomes 
scale free (Appendix |A| . 

While Equation ([3]) defines the amplitude of the veloc- 
ity on spatial scales /c"^, we need an actual velocity field 
to evolve our particle positions. To obtain it, we associate 
each shell m with three (fcmn, 'Umn) vector pairs indexed 
by n. These vector pairs characterize contributions to 
the gas velocity directed along Vmm and varying along 
kmn (for sheU m, |fc| = km)- We impose k^n -L *mn, 
so that the flow is incompressible. Further, the k^i and 
Vmi vectors for pair i of shell m are approximately per- 
pendicular to their counterparts in the other two pairs of 
that shell. This formulation (3 quasi-perpendicular vec- 
tors for both fc and i) for each shell in fc-space) allows us 
to span both fc-space and velocity space in each shell sep- 
arately while simultaneously avoiding strongly preferred 
directions. The projection of the velocity Vm associated 
with shell m onto the three vector pairs is done through 
the introduction of a unit vector a,„. The three com- 
ponents of am define the relative excitation of the three 
possible vector pairs. 
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The gas velocity is then given by 

2 ^ Ojm,n{i)Vm'^mn COs[kjnn ' ~t" <^mn(^)]; (4) 

where the outer sum is over shehs in /c-space, the inner 
sum is over the trio of {k,v) vector pairs, and 0,„„ is a 
phase shift included to avoid overlap at a; = 0. Note that 
amn is a scalar (the n-th component of am, itself a unit 
vector) while Vmn is the nth unit vector associated with 
the mth shell. Time dependence of the flow is allowed 
in three ways. First, the projection onto the k and i) 
vectors can be quenched (a™ and (j)mn held constant); 
these runs are named with a "Q" . Second, we can vary 
the unit vectors am for each shell through a random walk 
on timescales Tm, with (pmn held constant. This is our 
default method, with runs denoted by a "B" . Third, we 
allow the various phases (j)mn to vary as random walks 
while the am vectors are held constant in time; these 
runs are named with a "P". Note that the GOY model 
inherently includes phase variation by evolving a complex 
value for Vm (see Section l^.l.ip . 

We will use the subscripts Is (large scale) for rrimin 
and ss (small scale) for mmax- This approac h of using 
an art ificial velocity field is similar to that of iBec et al.l 
(j2005() . Our prescription is closer to that used in previ- 
ous analytical studies. It does however limit us to only 
considering turbulence in the inertial range as we have 
no means of tracking behavior and the energy injection 
or viscous dissipation scales. 

Our box size is L — 2^ and our boundary conditions 
are periodic, so the minimal k is k = 1 (m = 0). The 
largest k we will consider is fc = 256 (m = 8). In practical 
terms, assuming an energy-carrying scale at fc = 4 and a 
dissipative scale at a quarter of the Nyquist wave num- 
ber one would require approximately 8000'^ grid points 
to directly simulate the turbulence for our largest inertial 
ra nge (as compared to the very high resolution of 2048'^ 
in lBec et a l."2010al) . The benef it here can be seen in Fig- 
ure 1 oflCarballido et al.l (j2010t ) , where the inertial range 
is, perhaps, 10. Moreover, in their Figure 3, the differ- 
ence between the solid lines (assumed large Kolmogorov 
inertial range) and dotted lines (actual evolved turbu- 
lence) shows how large the effect on the particle motions 
of the limited inertial range actually is. 

2.1.1. GOY model 

The Gledzer, Ohkitani and Yamada (GOY) model 
(jOhkitani &: Yamadall 19891) solves an evolution equation 
for the energy in differing shells under the assumption 
that only shells adjacent in fc-space interact (to gener- 
ate a cascade rather than non-local effects): shells n, 
n + 1 and n + 2 can interact. We note that we solve the 
following system for a larger range of shells than are in- 
cluded in the synthetic turbulence that is applied to the 
particles: the choice m = (fc = 1) is associated with 
n — 4 in what follows. We follow N = 22 shells, using a 
slaved Ada ms-Bashforth scheme (jPisarenko et al.l 119931 : 
iMitra &: ba ndit 2004). Under this assumption we can 
write 

+ ^^n^ "^n = i {anV„+lVn+2+ (5) 
&„Z;„_lU„+l -I- CmVn-lVn-2)* + fn, 
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Fig. 1. — Velocity spectrum and time series for the GOY model. 
In the upper panel, we see that the power law of the GOY model is 
steeper than a Kolmogorov —1/3 law, closer to —0.4. In the lower 
panel, we show the velocity time series for the fc = 8 shell, with 
time normalized to that shell's own turn-over time. 

where w„ is the complex velocity associated with shell n, 
z/ is a viscosity and fn — for all n ^ 1 is the forcing 
term. The choice of constraints of energy and helicity 
conservation lead to 

a„ = kn, 6,1 = -fc„_i/2, c„ = -fc„_2/2, (6) 

with boundary conditions on the largest and smallest 
pairs of shells (n = 0, 1, iV - 1, A^) of 

&1 = 6^ = Ci = C2 = UN-l = flTv = 0. (7) 

The GOY shell mod el has a 3-cycle with shell index 
(jKadanoff et al.lll995[ ). which we filter using 

\vf\ = \lm{Vn+2Vn+lVn " W„+i /4) | (8) 

This model inherently provides time variation for the 
phase through the complex nature of u„. We combine 
this with Equation ^ by replacing Vm with \Vm\, set- 
ting (/)„i„ = arg(i;„i) and allowing a„i to vary as previ- 
ously described. 

In Figure [T] we show the velocity spectrum and a sam- 
ple shell velocity time series. We define a characteristic 
shell velocity Um through the time average of \um\ and 
Tm = ^/umkm, but this is not as perfectly defined as in 
the fully imposed Kolmogorov case. The GOY model 
velocity spectrum is also slightly steeper than a Kol- 
mogorov one, a consequence of its intermittency. 

2.2. Particles 

The motion of inertial particles in a fluid is determined 
by the friction between the particles and the fluid. The 
resulting drag force entrains the particles along the fluid 
motion. Particles with a finite friction time Tp (sub- 
script p referring to particles) are referred to as "iner- 
tial" and their motion deviates from that of the gas as 
long as the particles are not neutrally buoyant (ie. as 
long as Pparticie 7^ Pgas)- In this work we assume that 
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Fig. 2. — Log scale histograms of the total distance d traveled (in 
box units, box length 2tt) for 1000 particles as a function of time 
for a run with my^^^ = 2, mmax = 5. Overplotted is cfcj" (t/Tis)^^'^ 
where c = 2/3. 

Pparticic ^ Pgas which extremely well satisfied for proto- 
planetary disks. This allows us to neglect pressure forces 
on the particles. We initialize our particles with a fric- 
tional stopping time Tp, so that the equation for a parti- 
cle's velocity Up is determined by 



dup{t) _ Up{t) - V{xp,t) 
dt T„ 



(9) 



where V{xp,t) is the gas velocity at the particle's posi- 
tion Xp. It is this deviation of the particle velocity from 
that of the gas that allows particles to collide even in 
incompressible flows. 

When we evolve our particles in the velocity field given 
by Equation Q, we find a particle velocity dispersion 



a+Tp/tu)^/^'' 



(10) 



that matches the results expected since IVolk et ahl 
(jlQSOl ). Further, in Figure [5] we show log scaled his- 
tograms of the total distance d traveled in box length 
units (box side=27r, bin size=0.5) for a population of 
1000 particles in our base setup (see Section [S]). The 
time axis is scaled to the turnover time of the largest 
scale turbulence (time between samples 0.2tis). We over- 
plot 0.67fc;~^(t/r;s)^/^, which has the form and approx- 
imate scale of a random walk on the length and time 
scales of the largest scale turbulence. We can see that 
the particles' displacements are well fit by the random 
walk expected to be generated by the largest scale tur- 
bulent motion. 

In this work we consider only particles with identical 
stopping time Tp set to the turnover time of the shell nip 
with mmin < < TOmax- We wiU therefore refer to 
krrij, and Vmp as kp and Vp. In the limit of an infinite 
inertial range, this scale is the only scale, and so most 
of our values will be non-dimensionalized with respect to 
it. However, we can investigate the dependence of the 
collisional velocities on both smaller and larger eddies 



by turning on and off eddies of differing scales, i.e., by 
changing m,nin and mmax- 

We are primarily interested in dust collisions in proto- 
planetary disks, which allows some simplifications. The 
large size of the Kolmogorov scale turbulence in proto- 
planetary disks (at a minimum kilometer scale), along 
with the small size of the dust grains (quite sub-meter 
for grains embedded well within the inertial range) means 
that when considering the collisions of protoplanetary 
dust grains we must treat our dust-grains as point par- 
ticles, taking limits as the particle separation tends to 
zero. This is a distinguishing feature compared to at- 
mosp heric turbulence, whose Kolmogorov s cale might be 
1mm (|Shawll2003t iXu fc Bodenschatzll2008l) . comparable 
to rain drops or particulates. Finally, the actual collision 
rate of dust grains in protoplanetary disks is expected to 
be modest due to their dust grain number density: grains 
will experience multiple friction times between collisions. 
Accordingly, given our focus on protoplanetary disks we 
do not treat collisional scattering of the grains. 

It is traditional to non-dimensionalize the particles' 
stopping time to the Stokes number by scaling it to 
the turbulent turnover time at the dissipation scale (or 
occasionally to the time at the driving scale). Unfor- 
tunately neither definition fits our setup well, but one 
can define St ^ (kss/kp)^^-^ as a conceptual equiva- 
lent to scaling Tp to the dissipation scale since k^g is 
the smallest scale included in our simulation. Similarly, 
St ^ (kis/kpY^^ is a conceptual equivalent to scaling 
Tp to the turbulent time at the driving scale. However, 
in Appendix |^ we explain why we do not consider ei- 
ther approach to be conceptually important as we ex- 
plicitly hope to achieve large enough inertial ranges that 
the results are independent of St (either definition). We 
do define our effective Stokes number according to the 
first definition above, St' = (kss/kp)"^/^ , for use when we 
need to scale diagnostics to the smallest included scale. 
While it is not a perfect apples-to-apples comparison, 
this value sho uld be used when comparing wi th most 
other work fe.glBec et al.ll2010al : IPan et al.l[20TTI : but not 
lOrmel fc Cuzzill2007l where the other definition is called 
for). 

2.3. Run details 

In Table [T] we collect the details of the simulations we 
perform. The total inertial range is given by 2™"''''"™"'°, 
the effective Stokes number by St' = 22/3(™max-mp)^ ^nd 
the range between the largest scale turbulence and the 
particles is given by 2™p~™""". In all cases the time- 
step used to advance the particles was Clr^g. For runs 
using the GOY model, the GOY model is advanced for 
3 X 10^ time steps before the introduction of particles, to 
allow the turbulence to find its statistical steady state. 
The GOY equations are run with their own time-step 
of approximately 10"'* in code units (tweaked to fit an 
integer number of times within a particle's time-step); 
this shorter time-step is for the numerical stability of the 
algorithm. The particles are initialized with u = at 
random positions. 

We generate the k^n and Vmn vectors of Equation (jH) 
by randomly selecting kmi from the set of vectors with 
approximately appropriate length (the error decreases as 
the target |fc| increases) that fit in the periodic box. The 
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vector Vmi is generated as a random unit vector perpen- 
dicular to kjni- The vector km2 is then selected from 
the set of vectors of appropriate length that fit in the 
box and lie within 30 degrees of Vmi- Next we choose 
'Vm2 = kyni X kra2l\kmi X fcm2 1 , and repeat the process 
for fem3 and Vms- If no vector km2 or fcm3 can be found 
that satisfy the constraints, the process is restarted. 

3. ANALYSIS: BEHAVIOR 

We begin by making a full analysis of the run "Base" 
(Table [H smallest inertial range) to extract the behavior 
of the system, both with respect to clustering and colli- 
sion speeds, and to find a fit formula that can be applied 
in simulations of particle size evolution. In Section 2] 
we will study how both increasing the inertial range and 
changing our turbulence model affects the behavior. 

We will obtain our particle clustering and collisional 
velocity data by analyzing during run time snap-shots 
taken every full turbulent turnover tis for the largest 
eddy in the system; these snap-shots include every par- 
ticle's position and velocity. The particle positions are 
mapped onto a coarse grid, and every particle pair within 
a critical separation (at most 0.2 of the grid scale) is 
found. We bin our particle pairs simultaneously linearly 
in separation R, using either 50 or 100 bins depending 
on the run, and linearly in relative velocity u. We con- 
sider full spheres, so for every R we consider every par- 
ticle pair with \xi — X2\ < R- A bin listed with ve- 
locity Ubin contains particle pairs with relative velocity 
Ubin — Am < m < Mbin where Am = Umax/n uj ^max IS 
the largest considered particle pair relative velocity and 
n„ = 2000 is the number of velocity bins. For run Base, 
Au = 4 X 10""* in code units . Particle pairs with relative 
velocities larger than Wbin are included in that bin, so it 
is discarded. We will refer to the total number of pairs 
in bin (i?, u) of snapshot t as N(R, u, t), with dropped u 
implying summation over all velocity bins and dropped 
t implying averaging over snapshots. 

We took 529 snap-shots from run Base, but as dis- 
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for 


most of the analysis we drop the first 
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Runs 




Run 




rrip 


mmax Notes 








Basic Turbulence (Section [2TTJ 




B-A 





4 


6 Smallest Tp/rj^ 




B-B 





3 


5 




B-LI2 





3 


7 Largest inertial range kr/kg 


= 128 


B-C 


1 


3 


5 




B-LIl 


1 


3 


7 Large inertial range fey/fco 


= 64 


Base 


2 


3 


5 Subject of Section [3] 




B-D 


2 


3 


6 




B-E 


2 


3 


7 




B-F 


2 


3 


8 Largest St' = 10 






Quenched Spatial Projection (Section [JTlJ 




Q-A 


2 


3 


5 




Q-LI 





3 


7 Largest inertial range kr/kg 


= 128 






Varying Phase (Section 12. Ill 




P-A 


2 


3 


5 




P-LI 





3 


7 Largest inertial range kr/ko 


= 128 






GOY Turbulence (Section|2.1.1|l 




G-A 


1 


3 


6 




G-B 


2 


3 


5 




G-C 


2 


3 


6 




G-LI 





3 


7 Largest inertial range kr/ko 


= 128 



Rk =0.02 
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Fig. 3. — Particle pairs N{R, t) (unnormalized count) as a func- 
tion of time for run Base. Top panel: counts for a collisional sphere 
Rkp = 0.02. Middle panel: early time blow-up of top panel. Bot- 
tom panel: Clustering C(R), with a power-law fit in red/dashed. 

twenty, accordingly, N{R, u) and N{R) are averaged over 
509 snap-shots. 

3.1. Clustering 

In the top two panels of Figure |3] we show N{R, t) as a 
function of time for a representative maximal separation 
R. We further define the clustering C{R) as the ratio of 
N{R) to the expected number of particle pairs 

NiR)^-n,.^^^y (11) 

C{R) = N{R)/N{R), (12) 

assuming a spatially homogenous particle distribution, 
where Up is the number of particles and V is the volume 
of the box. A value of C{R) > 1 implies that parti- 
cles have been concentrated on length scales R, while 
C{R) < 1 implies some form of segregation or effective 
re pelling. The clust ering C{R) is the same as the g{St, r) 
of iPan et al.l (|2011[ ). modulo the differences between St' 
and their Stokes number. Their length scale rj might 
be understood as our kj^ , although their turbulence at 
that scale is explicitly affected by dissipation, and so the 
energy is no longer following a turbulent cascade. 

The particles are initialized with zero velocity and ran- 
dom initial position, and in the top and middle panels of 
Figure [3] we see that the number of particle pairs takes 
more than ten largest-scale turn-overs to stabilize (more 
than 16 particle friction times), and still shows signifi- 
cant fluctuations. This long stabilization time is an im- 
portant observation which helps explain results later in 
this paper: there is significant structure in the parti- 
cle distribution that takes time to develop, and studies 
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Fig. 4. — Average collisional speed between particles as a func- 
tion of maximal separation for run Base. Note the small scale of 
the horizontal axis. 

of particle collisions need to consider a significant time 
interval. Unfortunately, including a history of ten full 
largest-scale turnovers is beyond current or foreseeable 
analytical analyses. Clearly particle positions are quite 
correlated with one another. In what analysis that fol- 
lows, we discard the first twenty snap-shots to allow the 
system to stabilize. 

In the bottom panel of Figure [3] we show C{R) which 
shows a clear power-law dependence on the maximal sep- 
aration R. We will denote by fi the power- law C{R) oc 
and by i?i the intercept where the power-law fit 
gives 1. 

3.2. Average collision energy 

A first calculation of the rms-collisional velocity is 
shown in Figure |4] where we plot 

(^^N{R,u)v:- 

Uc[R) = 



^nN{R, u)u 



(13) 



The factor of u in the denominator is needed to convert 
from particle pairs to particle collision rates: because we 
are taking well-separated snapshots we need to weight 
high velocity pairs more heavily. A naive expectation 
would be that the collisional velocity decreases with de- 
creasing i?, reaching a finite plateau at small separation, 
such as seen in Bee et al. (2010b). This finite plateau is 
made possible by the compressible nature of the particle 
"fluid" , even if the turbulent gas itself is incompressible. 
In Figure S] we are already in this plateau, but see an 
increasing collisional velocity with decreasing R (except 
for the smallest R, where our pair-counts are too low for 
reliability). 

This perhaps surprising result does not flow from poor 
statistics, but rather is a robust and understandable con- 
sequence of the clustering seen in the bottom panel of 
Figure [31 The clustering seen there increases strongly 
with decreasing R, implying a relatively numerous pop- 
ulation of highly correlated particles. As R decreases, 
the relative particle-particle velocities in this highly cor- 
related population must decrease for the pairs to stay 
within the cluster radius: the smaller R, the lower the 
relative velocities of the correlated population must be. 
We might guess that the characteristic relative velocity of 
this population is linear in R, giving a constant crossing- 
time. Since the power-law fit in Figure [3] is weaker than 
i?^^, the contribution of this population to the denomi- 
nator of Equation decreases with decreasing R. The 
velocity plateau at larger separations implies that the 
highly correlated population has already ceased to con- 
tribute meaningfully to the numerator due to the ex- 
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Fig. 5. — Top and middle panels: particle relative velocity proba- 
bility distributions for three different R in run Base. Black: Rkp = 
0.02, blue/dashed: Rkp = 0.012, red/dash-dotted: Rkp = 0.008; 
see also the top panels of Figure [6] Bottom panel: the velocity of 
the peak as a function of R along with a linear fit in green. The 
vertical lines show the position of the three cases of the top two 
panels. 



tra factor of m^. Even at the largest separations R we 
consider in Figure IH the relative velocity of the highly 
correlated population is low enough that this population 
contributes only negligibly to the total collisional power. 
In effect, the highly correlated population is diluting the 
collisional energy when an average is taken, and this di- 
lution is decreased by considering smaller R. 

Since we are considering point particles, we must con- 
sider the limit i? — > 0, and in that limit, the highly corre- 
lated population appears to play a dominant role in the 
pair counts (because /i > 0), but contributes no collisions 
(since /i < 1). We accordingly define this coUisionless 
population as "cold" (low relative velocities). Unfortu- 
nately, numerical resource constraints prevent us from 
considering values of R small enough that the cold pop- 
ulation has ceased to contribute meaningfully to straight- 
forward calculations of the collision rate (the denomina- 
tor of Equation [T3| . Note that this is different f rom the 
expectations of, for example, iPan et aD ()2011| ). where 
turbulent clustering was hypothesized to lead to more 
common collisions by increase the number of nearby po- 
tential collisional partners. As we will see, however, there 
is a significant chance that some level of clustering will 
persist, albeit not one that scales inversely with R. 

The implications of this interpretation of Figure |4] are 
seen in Figure El where 



N'{R,u) = N{R,u)/Au 



(14) 
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Fig. 6. — Top panels: fits using Equation l|16|l for run Base. 
Black/dash-triple-dotted: data, red/solid: total fit, green/dashed: 
fit component a, blue/dash-dotted: fit component b. Only the 
region to the right of the vertical line u = Uc is fitted. Left panels: 
Ua and ui,, along with their maxima in the non-contaminated region 
to the left of the vertical line R = Re- Right panes: Ca and Ct, 
along with their means in the non-contaminated region to the left 
of the vertical line R. = i?,c . 

is the density of pairs within separation R in velocity 
space (to convert between binned data and a smooth 
distribution). The probability distribution is seen to be 
strongly affected by a low velocity peak which moves to 
smaller velocit y linearly with J?, to a limit of as i? 
(compare with iCarballido et al.ll20fol Figure 5, which is 
lin-log instead of log-lin). We define 



UJ = Ucold/R, 



(15) 



as the measure of the dependence of the velocity of the 
cold population on separation, where Ucoid is the modal 
pair velocity. 

This collisionless behavior of the cold population does 
not mean that we expect no collisions. As can be seen 
in the top panel of Figure [SJ N{R, u) has a significant 
high velocity tail, which can contribute "real" collisions. 
A process to excise the cold, collisionless population and 
fit the remainder is discussed below. 

3.3. Fitting 

We have failed to successfully fit the cold, clus- 
tered population with a tractable distribution such as 
a Weibull distribution; however the tail of the pair- 
probability distribution in velocity space is well fit by 
a simple exponential, and at more modest velocities the 



addition of a second exponential improves matters. Ac- 
cordingly, we fit the particle-pair probability distribution 
by 



N'{R,u > u,)/N[R) = ^e-"/"" + ^g""/"' 

Ua Ub 

Uc = O.lWr, 



(16) 
(17) 

but only for u > Uc, where Uc is a cut-off designed to ex- 
clude the cold population. The fit formula Equation (|16|) 
identifies two separate populations a and b with Ua < ut- 
While we do not have an analytical theory for the form of 
Equation (|16p , the quality of the high velocity component 
b can be seen from the overlap between the blue/dash- 
dotted curves (the higher velocity component) and the 
black/dash-triple-dotted curve (data) in Figure HI The 
fit for component a to be less well constrained, but it 
gives a noticeable different at lower relative velocities, 
seen in the difference between the blue/dash-dotted curve 
(only component b), and the red/solid curve (both com- 
ponents, overlays the data outside of Uc). We refer to 
component a, handling the lower relative velocities, as 
the "warm" population, and the fit component b, han- 
dling the higher relative velocities, as the "hot" tail. 

The parameters Ca and Ct then set the clustering of 
the effective densities of the warm and hot populations 
compared to the total population were it homogeneously 
distributed (see Equation [12]). Finally, we define a con- 
tamination radius 



Rc ^ 0.2uao/uj 



(18) 



outside of which the cold population is considered to con- 
taminate the results, i.e., where Ucoid from Equation (jlSp 
is comparable to Uq. The factor of 0.2 is chosen ad-hoc, 
but is inspired by the middle panels of Figure [6l Uao 
is an approximate value for Ua- The fit parameters in 
Equation (fT6|) are functions of R, but as seen in Figure[6] 
meaningful limits i? — t- can be taken. The reported val- 
ues of Ca, Cb, are averaged over the non-contaminated 
region, while the reported values of Ua, and Ub are the 
maximum value inside the non-contaminated region. 

3.4. Weights and errors 

Our reported fit coefficients Ua,Ub,Ca and Cf, depend 
on the weighting scheme, the fitting region, and Uao- We 
weight with 



W = [l/N{R,u)] X 



1 



2(it — 0.7um)('" — Uc) 



(0.7uM 



, (19) 



where um is an upper velocity cut-off imposed when the 
number of pairs/bin approaches unity. This weighting 
scheme weights the low velocity (warm) and high velocity 
(hot) regions more highly, allowing the two populations 
to be distinguished by the fitting algorithm. The values 
of N'{R,u) in Figure [6] have a minimum ^ 10. As we 
are averaging over 509 snap-shots, with Au = 4 x 10^^, 
that corresponds to ~ 2 pairs per bin which is a conse- 
quence of our large number of velocity bins, 2000. To 
avoid zeroes we have smoothed over 20 bins in velocity 
space. The effect of sampling noise can be gauged from 
the middle and bottom panels of Figure [S] as the limit 
i? — > is taken with a corresponding reduction in pair 
count, especially at high relative velocity. As the fitting 
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procedure is somewhat qualitative, we will be more in- 
terested in the scale of the results, and any trends with 
varying TOmin and TOmax, rather than in precision. As we 
will discuss in Section [5l our collisional velocities are fac- 
tors of order five slower than previous work, well outside 
any systematic or statistical errors. 

4. ANALYSIS: INERTIAL RANGE, TURBULENCE MODELS 
4.1. Clustering 

In Table [2] we present the basic clustering data for our 
turbulence models, across an array of inertial ranges and 
our turbulence models. Parameters TOmin, Wp, rrimax and 
St' define the included inertial range, as well as an ef- 
fective Stokes number. The diagnostics fj, and Ri are 
the power-law and intercept of the power-law fit to C (R) 
(Section l3.1|) while uj defines the relative velocities of the 
cold population (Equation [T5|) and Rc defines the con- 
tamination radius outside of which the collisionless cold 
population becomes indistinguishable from the collisional 
warm and hot populations. 

There appears to be a modest decrease of the cluster- 
ing parameter /x with the number of shells with turnover 
times shorter than the particle drag time Tp, as mea- 
sured by St' . However, it also appears to increase with 
the number of shells included with turnover times longer 
than Tp. Comparing runs B-E, B-F, B-LIl, and B-LI2, we 
feel that we have reached a large enough inertial range to 
estimate that particles well embedded in turbulence will 
experience clustering with /x ^ 0.4. This result matches 
the behavior with other turbulence models, although run 
P-A is a clear outlier. 

Another interesting column is ujTp. This value scales 
slightly slower than St' , and oj"^ is the crossing time of a 
cold cluster. Accordingly, the result is that the crossing 
time of the cold cluster scales slightly slower than the 
turnover time of the smallest included eddy, supporting 
the extension of the hypothesis that the cold velocity is 
linear with separation (Eq. I15|) to an infinite cascade 
without meaningful smallest scale. This is further born 
out by the results for the other turbulence models, partic- 
ularly the GOY model, which see weak or no dependence 
of oj on St'. 

A problematic column is Rckp. This imposes the min- 



TABLE 2 
Clustering diagnostics 



Run 


mini 


in rup 




St' 




Rikp 




Ftckp 


B-A 





4 


6 


2.5 


0.65 


1.6 


0.29 


0.04 


B-B 





3 


5 


2.5 


0.57 


1.9 


0.27 


0.04 


B-LI2 





3 


7 


6.3 


0.39 


1.6 


0.55 


0.02 


B-C 


1 


3 


5 


2.5 


0.55 


1.7 


0.23 


0.04 


B-LIl 


1 


3 


7 


6.3 


0.38 


1.8 


0.55 


0.02 


Base 


2 


3 


5 


2.5 


0.49 


1.6 


0.29 


0.03 


B-D 


2 


3 


6 


4 


0.45 


1.5 


0.39 


0.02 


B-E 


2 


3 


7 


6.3 


0.38 


1.2 


0.52 


0.02 


B-F 


2 


3 


8 


10 


0.38 


1.2 


0.76 


0.02 


Q-A 


2 


3 


5 


2.5 


0.34 


5.0 


0.51 


0.03 


Q-LI 





3 


7 


6.3 


0.44 


13 


0.50 


0.03 


P-A 


2 


3 


5 


2.5 


0.81 


0.8 


0.21 


0.03 


P-LI 





3 


7 


6.3 


0.40 


2.2 


0.39 


0.03 


G-A 


1 


3 


6 


3.7 


0.36 


1.9 


0.14 


0.02 


G-B 


2 


3 


5 


2.4 


0.40 


2.0 


0.16 


0.03 


G-C 


2 


3 


6 


3.7 


0.33 


2.7 


0.14 


0.03 


G-LI 





3 


7 


5.9 


0.38 


1.8 


0.20 


0.02 



imal spatial resolution required to make a meaningful 
analysis of dust collisions by discarding the cold popula- 
tion. In a 27r'^ box with kp = 8 (assuming an energy car- 
rying scale of fc = 4), a grid scale with Ax — Rc requires 
a resolution on the order of 1700'^ even for our smallest 
inertial range. This smal lest range is within reach, albeit 
barely (|Bec et al.ll20T0al|O ). 

4.2. Collisions 

In Table[3]we collect fit data for our simulations, which 
should only be used for intermediate and higher relative 
velocities {u > Uc = O.lvp, see Equations [TBI and [T7)) . Pa- 
rameters mjiiin, TUp, minax and St' define the included in- 
ertial range, as well as an effective Stokes number, while 
the diagnostics Ca,Cb,Ua and are the coefficients in 
our fit formula Equation p6|l . Coefficients Ca and Cb 
describe the effective target number density seen by the 
warm and hot populations respectively, and coefficients 
Ua and Ub represent the velocity scale of the two popula- 
tions. 

The Ca and C'b columns indicate that the warm and 
hot populations, in the velocity regime where the fit ap- 
plies, have individual effective densities comparable to 
the density that the total dust population would have 
were it evenly distributed. This does not strictly imply 
that there is clustering of the collisional population, be- 
cause many of the collisions are likely at low, unfitted 
velocities. The results for Ca and Va are the least useful, 
as they are most sensitive to the cold population, and 
Ua < Uc- The results for Ub are nicely only weakly sen- 
sitive to the set of included shells, while Cb appears to 
require a significant St' ~ 10 to reach a limiting value. 

Most of our diagnostics, such as Ua,Ub, Ca and i?i ap- 
pear to be relatively insensitive to the set of included 
shells, i.e., (rTimin — mmax) and St' . It should nonetheless 
be noted that since the diagnostic Ub is crucial to the 
high velocity collisions which might result in fragmen- 
tation, its weak sensitivity to mmin is still important. 
Achieving converged values for C(,, an important diag- 
nostic, certainly requires St' > 6, while the behavior of 
Ca in the limit of both larger and smaller eddies is not yet 
clear (bottom two lines of Table [3]). Accordingly, ranges 
of kp/kis — 8 and k^s/kp ~ 16 appear to be reasonable 



TABLE 3 
Collisional diagnostics 



Run 


^min 


rrip 




St' 


Ca 


Cb 


Ua/Vp 


Ub/Vp 


B-A 





4 


6 


2.5 


1.1 


0.46 


0.070 


0.21 


B-B 





3 


5 


2.5 


1.2 


0.57 


0.056 


0.20 


B-LI2 





3 


7 


6.3 


1.7 


0.92 


0.069 


0.19 


B-C 


1 


3 


5 


2.5 


0.9 


0.42 


0.054 


0.18 


B-LIl 


1 


3 


7 


6.3 


1.5 


0.89 


0.069 


0.17 


Base 


2 


3 


5 


2.5 


1.0 


0.62 


0.046 


0.15 


B-D 


2 


3 


6 


4 


1.1 


0.83 


0.054 


0.17 


B-E 


2 


3 


7 


6.3 


1.1 


1.0 


0.057 


0.16 


B-F 


2 


3 


8 


10 


1.2 


1.1 


0.064 


0.17 


Q-A 


2 


3 


5 


2.5 


1.4 


1.4 


0.077 


0.48 


Q-LI 





3 


7 


6.3 


2.8 


2.0 


0.13 


0.71 


P-A 


2 


3 


5 


2.5 


0.5 


0.25 


0.036 


0.098 


P-LI 





3 


7 


6.3 


1.5 


0.37 


0.050 


0.134 


G-A 


1 


3 


6 


3.7 


0.7 


0.31 


0.070 


0.22 


G-B 


2 


3 


5 


2.4 


0.7 


0.43 


0.077 


0.32 


G-C 


2 


3 


6 


3.7 


0.7 


0.31 


0.068 


0.20 


G-LI 





3 


7 


5.9 


0.8 


0.37 


0.076 


0.24 
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minimums. Unfortunately, it appears that extreme res- 
olution is required not merely to include an adequate 
inertial range, but also to capture the details of particle 
clustering. 

Comparing the results across different turbulence mod- 
els, a first observation is that they are well fit by Equa- 
tion (ITB)) . as seen in Figure [71 albeit with differing pa- 
rameters. The fully quenched model (runs Q-A and Q- 
LI) produces both significantly more concentrated colli- 
sions (the C coefficients), and significantly higher col- 
lisional velocities. This is not particularly surprising 
as there is nothing to destroy long time correlations. 
Putting the time variation of the spatial projection into 
the phase (runs P-A and P-LI) however results in a sig- 
nificant drop in the collisional velocities. Finally, the 
GOY model appears to fit with very similar velocity pa- 
rameters of our base model although the clustering is 
weaker. This is somewhat surprising as the turbulence 
is spiky, which might be expected to result in a stronger 
high velocity tail. On the other hand, it inherently in- 
cludes phase variation. A further unexpected detail of 
the GOY model is that the smallest range of included 
shells (TOmin = 2, TOmax = 5) showed the strongest high 
velocity tail {ui, — 0.32wp). This may be due to the fact 
that the GOY scheme has an intermittent energy cas- 
cade, and the generation of sub- vortices from larger vor- 
tices (a consequence of localized energy spikes in fc-space, 
through the coupling terms of Equation [5]) , if appropri- 
ately followed by including an adequate inertial range to 
track the sub- vortices, will decrease particle correlations. 

Finally, we suggest that Equation be used with 
parameters of the order of Ca ~ 1-5, Cf, ~ 1, ~ 0.07vp 
and Ub — 0.2vp, applying only when relative velocities 
exceed the cut-off velocity Uc = O.lvp, i.e.. 



Rk =0.014 



Rk =0.020 



N{R,u> O.lvp) 

Wr) 



(20) 



If one prefers the GOY model, the concentration param- 
eters should be changed to Ca ^ 0.8 and Cb ~ 0.4. 

5. COMPARISON WITH PREVIOUS WORK 

5.1. Clustering 

Turbulent clustering has been a topic of interest in 
the protoplanetary community, and has been invoked, 
for example, in the creation of chondritic meteorites 
(jCuzzi et al.l [^01). Our clustering diagnostics fol- 
low those of Pan ct al. (2 0Tll). rat h er th an those of 
iFessler et all ([1994) : Hogan fc Cuzzil (|200l . The pre- 



vious work found evidence for strong clustering for par- 
ticles with modest St, i.e., particles reasonably well cou- 
pled t o turbulence at t he dissipation scale. The cluster- 
ing of IPan et al.l (poTTI ) was observed to decrease as their 
Stokes number became di fferent from 1. O ur results fit 
in magnitude with those of lPan et al.l ()2011|) , which scale 
up to ^ > 0.6 (see their Figure 6). However, our cluster- 
ing is significantly stronger at high St or St' . Given the 
differences in the fluid flows used (synthetic turbulence 
in our case, forced turbulence in iPan et al.i i20lTI) , the 
similarities in behavior and scale of the exponent are en- 
couraging. However, it isn't clear from the /x column in 
Table [21 that clustering will cease for large St' , as previ- 
ously seen. It should be noted that previous simulations 
had limited inertial ranges, so their high St particles are 
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Fig. 7. — As in Figure [6] but here wc show run G-LI. Averaged 
over 112 snap-shots, with Am = 4.64 X 10^* in code units. 

not fully embedded in the turbulence (see Appendix 1X1) . 

There is reason to expect that clustering will be more 
effective for particles with St = 1 (Cuzzi ct al. 200il; 
iHogan fc Cuz^l200l . However, we can imagine evolv- 
ing a system of two species of particles with Sti = 1 
and St2 <C 1. Once a statistical steady-state is achieved, 
with the Sti particles more strongly clustered than the 
St2 particles, we decrease the fluid viscosity significantly 
so that in the new fiow, Sti ^ 1 and St2 = 1 because 
the new dissipative scale is smaller. At that point, even 
if the St2 particles are more strongly clustered than the 
Sti particles, it isn't clear why the Sti particles would be 
less clustered than originally. Runs Base through B-F in 
Table [21 match this thought experiment, and do show a 
drop of fi from 0.49 to 0.38 which appears to be a minimal 
value. This weak decrease in /i may result from particle 
clusters being disrupted by the additional, smaller-scale, 
e ddies. 

IPan et all (|2011f) hypothesized that clustering, espe- 
cially if it exists for St > 1, might greatly increase the 
collision rates between particles by generating regions of 
enhanced particle number density. Our results allow us 
to contradict the straightforward version of that hypoth- 
esis: the clustering exists, but the particle clusters are 
non-collisional. However, once particle number densities 
get high enough that the dust fluid density is compara- 
ble to the gas density, the dust drag has a signiflcant 
back-reaction on the gas. An example is the streaming 
instability. In protoplanetary disks, gas orbits in a sub- 
Keplerian fashion because of the outwards pointing pres- 
sure force. As a result, particles, which would naturally 
orbit in a Keplerian fashion, feel a headwind. Similarly 
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to drafting on highways or in bicycle races, clumps of 
particles can then form which are dense enough to back- 
react on the gas (iCuzzi et al.ll2001l: lYoudin fc GoodmanI 
120051: iJohansen et al.ll2007t iLewellen et alJl2008l )! 

5.2. Collisions 

Much of astrophysical turbul ent dust c ollisio n 
work follows th e app ro ach of IVolk et al.l (1980[ ): 
Markiewicz et all (I1991D : ICuzzi fc HoganI (|2003l) :~^e 
Ormel fc Cuzzil (1200711 for a recent. This approach de- 
fines two classes of eddies, Class I large-scale eddies and 
Class II small-scale eddies. The former have large veloci- 
ties and time scales, and so can transport dust grains sig- 
nificant distances: they dominate the turbulent diffusion 
of dust throughout a disk or atmosphere. On the other 
hand, they change slowly enough in both space (com- 
pared with dust stopping lengths) and time (compared 
with Tp) that nearby dust grains, which could collide, 
see nearly identical gas motion. As a result. Class I ed- 
dies can affect the coUisional behavior of dust grains only 
slightly. Class II eddies on the other hand vary rapidly 
compared to both the frictional stopping time of the dust 
grains and their stopping length. Accordingly, these ed- 
dies can affect even nearby dust grains differently, driving 
collisions. However, their short time and length-scales 
mean that they provide only weak large-scale transport. 

Dropping the contribution of Class I eddies for col- 
lisions between identical dust grains, lOrmel fc Cuzzil 
(|20?il find an rms collisional velocity of 



Extended 



Flat 



^Ormol 



(21) 



While this result is frequently quoted in terms of the 
Stokes number defined as St = Tpftx, the above ver- 
sion is more general (see Appendix An important 
assumption is that the relative motion of particles ap- 
proaches a finite limiting value as their separation goes 
to zero, which is possible for inertial particles (r^ ^ 0) 
whose motion must deviate from that of the gas. Such 
behavior i s seen in the direct numerical simulations of 
iBec et all (|2Q10bf ) and we believe we understand the de- 
viation from that assumption seen in our Figure |4l This 
analytical approach cannot however handle very long 
time-correlations between dust grains. 

The result in Equation (|2ip is different from our Equa- 
tion in two ways. First, it is a single number, which 
does not suffice to describe our results, with their two 
velocity scales of the warm and hot populations. Our 
results also have merely exponentially falling probabil- 
ity tails with u, a much broader distribution than, for 
example, a Maxwellian distribution. They should there- 
fore be treated as a probability distribution. Second, if 
one nevertheless uses Equation (fT3)) to extract a single 
rms-averaged collisional velocity, we find 



6 {CaUl 



1/2 



WOrmcl/4, (22) 



using the suggested parameters of Section 14.21 a signif- 
icant decrease in the characteristic collisional velocities 
and a much larger decrease in the turbulent collisional 
energies. We attribute this difference in predicted colli- 
sional velocities to the high level of correlation we have 




0.0 0.2 0.4 0.6 0.8 




Fig . 8. — Ad hoc fits for u < Uc- Left column: extending Equa- 
tion ||20|I to ii = 0, Riglit column: extending tiie constant valu e 
at -u = Mc to M = 0. Top panels: pair density, i.e., Equation I I20I I. 
Middle panels: solid: cumulative totals of the top panels; dashed: 
Uc its intercepts by the data. The right axis is self-normalized. 
Bottom panels: solid: self-normalized cumulative tot als o f the col- 
lision rate i.e., the first velocity moment of Equation II20I I: dashed: 
Uc its intercepts. 

identified, which not only creates the "cold" population, 
but appears to slow all the collisions. 

One can also compare Figur e ID and the botto r n pane l 
of Figure [5] with Figure 4 of iCarbalhdo eTall (|2010t) , 
which is further limited by the sub-gridscale gas veloc- 
ity interpolation. Clearly, very high spatial resolution 
is required to extract actual collisional velocities from 
numerical simulations, at least for particles with identi- 
cal stopping times. Our method provides large inertial 
ranges and high resolutions, whose necessity can be seen 
in their Figure 3, where the analytical result based on a 
full inertial range, the analytical result based on the ac- 
tual inertial range and the numerical result all disagree. 
The resolution to begin to see the multiple populations 
we have identified in direct numerical simulations of the 
full Navier-Stokes equation may be becoming availab le 
(for example the 2048^ simulation of IBec et alll2Q10allbl 
which is just adequate to include the inertial range as- 
sumed in run Base while resolving Rc). However the 
analysis has not been done in the same terms so it is not 
clear whether they would see the effects we predict, or, if 
not, it would be because they contradict our results, or 
do not have adequate resolution. 

6. DISCUSSION AND CONCLUSIONS 

We arrive at a picture of turbulence-induced parti- 
cle concentration and collisions which creates spatially 
small (with respect to turbulent scales), highly corre- 
lated "cold" clusters of particles. While these clus- 
ters are themselves non-coUisional, they periodically pass 
through each other, resulting in bursts of collisions. 
These clusters are channeled through the high-strain 
boundaries between turbulent eddies, and even separate 
clusters are strongly correlated. This correlation causes 
the turbulent collision speeds we find to be significantly 
slower than those predicted by analytical approaches, 
that cannot treat the long time-scale correlations. This 
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picture of particle-particle collisions occurring when clus- 
ters of particles pass through each other is supported 
by the long stabilization time seen in Figure [H where it 
takes many turbulent turnover times to reach a statistical 
steady state, implying non-trivial dust spatial structures. 

Our results apply most strongly for particles deeply 
embedded in a turbulent cascade. Our approach does 
not have the ability to handle the details of the forcing 
scale or the dissipation scale of the turbulence, both of 
which could have different behavior suc h as the bottle - 
neck effect at the dissipation scale (Dob ler et al.l [20031 ). 
Nevertheless, the qualitative similarity between our low 
inertial range runs (Basic, Q-A, P-A, G-A) with the large 
inertial range runs (B-LI2, Q-LI, P-LI, G-LI) suggests 
that this will only pose a difficulty if the statistics of the 
turbulence are indeed quite different at those two scales. 
A second difficulty is that our synthetic turbulence model 
cannot handle the dust's back-reaction on the gas, which 
can be significant if there is strong clustering. An im- 
plicit assumption of our work therefore is that the dust 
spatial density is less than that of the gas. 

We find that Equation ([T51) is a workable fit for the 
number of turbulence-induced particles pairs as a func- 
tion of coUisional velocity. This form has some interesting 
implications. A simple one is that it is broad: unlike the 
case of, for example, a Maxwellian distribution, enough 
of coUisional velocity space is populated that many col- 
lisional outcomes will occur in relevant quantities. Less 
obviously, we note that in the "high" velocity parameter 
regime, we can fit the coUisional probability distribution 
through the autocorrelation of exponential decays: 
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This implies that (faster) collisions are generated by par- 
ticles each contributing a velocity with respect to the col- 
lisional center (which is a feature of the fluid flow rather 
than the center of mass of the particle pair) with an ex- 
ponentially falling probability. Note that this is not the 
difference between the particle velocity and that of the 
gas, which is, naturally, correlated between the two parti- 
cles. Why the individual particle probability distribution 
is an exponential decay, and why the autocorrelation oc- 
curs in collision-space {P{u) oc ite^"/"'') rather than in 
the particle pair space {P{u) oc e^"/"") is unclear. 

Our formula is limited in its applicability to m > u^, 
and as R tends to 0, neglects almost all the particle pairs 
since they are part of the cold population (because /i > 0, 
see the height of the peak in the middle panel of Figure [SI 
which increases as R decreases). This does not appear 
to pose unsurpassable difficulties: in Figure [8| we show 
two alternate ad-hoc extensions for Equation ([20]) into 
the low u-regime. From the bottom panels we can see 
that it is unlikely that errors in the collision rate based 
on the extension chosen will be more than 25% of the 
total: low relative velocity dust grains, even if numerous, 
contribute few collisions due to their slow speeds. Errors 
in energy, the square-root of the 3rd velocity moment, 
will be significantly less yet, as long as there is not a 
hidden, highly clustered but still coUisional population 
that we have not been able to detect so far. 

Another intriguing result is the low scale of Ub — 
0.2vp. This coUisional velocity scale is significantly be- 



low the predicte c i rms coUisional velocity = 2i/2„ of 
lOrmel fc Cuzzil ()2007[ ). This might be explained by the 
vertical axes of the middle plots of Figure [Sj the limiting 
value on the left axes gives the effective concentration 
i.e., the number of pairs in the warm and hot popula- 
tions, divided by N{R). While the difficulties fitting the 
cold population make the result uncertain, it appears 
that clustering occurs in the coUisional population: par- 
ticles see more warm and hot coUisional partners than 
would be expected if all dust grains were homogeneously 
distributed in space. This clustering is different from 
that of the cold population in that the clustering is in- 
dependent of R for small separations (i.e., the effective 
^ = 0; an equivalent plot to the bottom panel of Figure [3] 
would be flat). For any clustering to occur, the collid- 
ing particles in the warm and hot populations must still 
be significantly c orrelated with one ano ther, an aspect 
which the work of lOrmel fc Cuzzl (j2007D could not fully 
capture. 

In conclusion, we see turbulence-induced dust colli- 
sions at significantly (a factor of 3-5) slower velocities 
than predicted by existing analytical theory. If we con- 
sider an approximate a-disk Minimum Mass Solar Neb- 
ula with a sound speed of Cg = lO^cm/s (temperature 
T ~ 300 K) at lAU, a = 10~^ and centimeter-scale 
grains with a stopping time Tpflx — 0.01 and mass 
m ^ 1 g, the turbulen t velocity at the larges t scale is 
yt = 3 X lO^cm/s ( Shakura fc SunvaevI 119731 : iHavashil 
|1981| ). Un der these condit i ons, t he predicted coUisional 
velocity of lOrmel fc Cuzzil (j2007[ ) is m = 5 x 10^ cm/s, 
well in to the destructive collision regime of lGiittler et al.l 
(|2OT0h . Fi gure 11. A five- fold reduction in the collision 
speed would however lead to bouncing. While this result 
implies frequent destructive collisions, we also predict 
large numbers of lower velocity collisions, which could 
lead to sticking. This would ease the difficulties in plan- 
etesimal formation associated with bouncing and frag- 
mentation (Zsom et al. 2010). The velocities we expect 
will still lead to significant fragmentation, but the in- 
clusion of a slower coUisional velocity probability distri- 
bution allows for the consideration of "lucky" particles 
that are in unusually low velocity collisions to grow large 
enough that they can survive. While the slower coUi- 
sional velocities reduces collision rates, we also see a lim- 
ited enhancement of the effective dust number density 
through clustering, which mitigates this collision rate re- 
duction. 

Finally, the hot population coUisional tail falls off only 
exponentially with velocity, rather than something closer 
to a Maxwellian distribution. This emphasizes the im- 
portance of using a coUisional velocity probability dis- 
tribution instead of a single characteristic coUisional ve- 
locity. For example, high velocity outlying events are 
expected to occur at non-negligible rates, and could con- 
tribute to a fragmentation cascade and dust reprocessing. 
This would occur even if the reduction in coUisional ve- 
locities results in few enough destructive collisions that 
the fragmentation barrier to dust growth is lifted. 

Our results are well suited to inclusion in a model of 
coUisional dust grain agglomerat i on in protoplanetary 
disks such as that of I Zsom et al.l (|2010f). The formula 
given by Equation (jl6p is simple enough for inclusion, 
while allowing a velocity probability distribution. This 



12 



enables the full use of experimental results about the crit- or fragment, 
leal velocities at which colliding particles stick, bounce 

APPENDIX 
STOKES NUMBERS 

Studies of turbulent particle transport generally use the Stokes number St = Tp/r^ to non-dimensionalize the 
stopping time, where r,, is the turbulent time scale associated with the viscous dissipation scale. Astrophysical 
studies of protoplanetary disks however often use St = Tpflx, where is the Keplerian rotation rate, since it is 
better constrained than and is believed to be a good estimate for the time scale associated with the largest scale 
turbulence. As such, the particle stopping time is often scaled to two completely different time scales, the largest 
and the smallest associated with the turbulence. Neither formulation is however clearly appropriate for studies of 
particle-particle relative motion and clustering. 

Stating that the relevant quantity is St = Tp/rrj implies that either the details of the dissipation process proper, or 
the lack of fluid energy at scales fe > fc^ plays a crucial role in the particle response to the turbulence. For a particle 
with St <1, such effects are expected, as particles couple most strongly with eddies with turnover times t ^ Tp and 
the dissipative cutoff means that many of those eddies are missing. However, it is unclear why a dependence on Tp/r,, 
would exist for particles with St ^ 1. Instead, for particles with rjs ^ Tp ^ r,,, we expect the behavior to be scale 
free: the largest scale of the turbulence (r/^) is too large to affect particle-particle relative motion while the dissipation 
scale is too small, so the turbulence does not set time or length scales for the particles. In this regime, instead of 
measuring the particle stopping time as a function of some turbulent t ime, one sh o uld in stead measure the turbulence 
by the particle stopping time. This approach is implicitly followed bv iVolk et al.l ()1980D . where the division of eddies 
into Classes I and II is done by measuring the turbulent turnover time relative to Tp. They do quote results as a 
function of St = TpJ7i<-, but that is only possible because of the assumption of a Kolmogorov cascade, which allows 

them to note that eddies with turnover times t — Tp have velocities Vp — y/Sivis- 

This poses difficulties in numerical simulations of particle relative motion in turbulence because the accessible ratio 
of Tis/Trf is modest at best. Accordingly, any apparent dependence on St = Tp/rrj is difficult to distinguish from a 
dependence on St = Tp/ris and even then would i mply little for th e case of a particle deeply embedded in a large 
cascade. For example, when considering the work of lPan et al.l ()2011[ ). it should be noted that their ability to track the 
clustering of particles deeply embedded i n a turbulent cas c ade (x p ^ r^s but also Tp <^ tis) is limited by their modest 
inertial range. Alternatively, Figure 3 of iCarballido et all ()2010D plots the difference between theoretical results for a 
large inertial range (solid) and the theoretical results applied to the numerically obtained inertial range (dotted) . Any 
extension of their results of decreased clustering for large St is suspect for a protoplanetary disk with a large inertial 
range. It does, however, certainly imply that clustering will decrease once Tp/r^s ~ Tpflx becomes large. 
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TABLE 4 

Variables, Parameters an Diagnostics 



NamG 


IN U LCb 


Section 


^min 
^max 


Shell numbering for largest scale shell 
Shell numbering for smallest scale shell 
Shell numbering for shells with k = kp 


Section O 


Is 
ss 

p 


Subscript for large-scale, equivalent to mmax 
Subscript for small-scale, equivalent to mmin 
Subscript for particle, equivalent to nip 


Section [2n 
Section [Ml 


V{x,t) 

Tla 

"J" 3 3 


Gas velocity 

Turbulent turnover time at the largest turbulent scale 
Turbulent turnover time at the smallest turbulent scales 


Equation Q 
Section EH] 


i^mn 
l^mn 
"^mn 
dmn 
4^mn 


Synthetic turbulence parameters 


Section O 


Up 
Tp 

hp 

Vp 

St' 


Particle velocity 
Particle stopping time 

Wavenumber of turbulence with turnover time t = Tp 
Gas velocity at fc = fcp 
Effective Stokes number 


Section [12] 


Ua 
Uf, 
Uc 
Ca 
Ci 


Coefficients for the particle-pair number and velocity 
probability distribution fit 


Equations II16II and II17II 






N(R,u,t) 

N{R,u) 

NlR,t) 

N{R) 

CiR) 

Ri 


Particle pair count within separation R, at velocity u in snap-shot t 
as N{R,u,t), averaged over snap-shots 
as N{R,u,t), integrated over velocity 

as N{R,u,t), integrated over velocity and averaged over snapshots 

Particle concentration 

Power-law exponent for C(R) 

y = 1 intercept for power-law fit of C(R) 


Section [XT] 


N'(R, u) 

LJ 

Rc 


Pair density in velocity space: N{R,u)/Au 
Cold population velocity measure 
Contamination radius 


Equation (14)1 
Equation II15II 
Equation II18II 
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